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ABSTRACT 


The next generation gravitational waves (GW) detectors are most sensitive to 
GW emitted by compact (neutron star/black hole) binary mergers. If one of those is 
a neutron star the merger will also emit electromagnetic radiation via three possible 
channels: Gamma-ray bursts and their (possibly orphan) afterglows (Eichler et al. 
1989), Li-Paczynski Macronovae (Li & Paczyhski 1998) and radio flares (Nakar & Pi- 
ran 2011). This accompanying electromagnetic radiation is vitally important in con¬ 
firming the GW detections (Kochanek & Piran 1993). It could also reveal a wealth 
of information regarding the merger and will open a window towards multi-messenger 
astronomy. Identifying and characterizing these counterparts is therefore of utmost 
importance. In this work we explore late time radio flares emitted by the dynami¬ 
cally ejected outflows. We build upon previous work and consider the effect of the 
outflow’s non-trivial geometry. Using an approximate method we estimate the radio 
light-curves for several ejected matter distributions obtained in numerical simulations. 
Our method provides an upper limit to the effect of non-sphericity. Together with the 
spherical estimates the resulting light curves bound the actual signal. We hnd that 
while non-spherical geometries can in principle lead to an enhanced emission, in most 
cases they result in an increase in the timescale compared with a corresponding spher¬ 
ical conhguration. This would weaken somewhat these signals and might decrease the 
detection prospects. 


1 INTRODUCTION 

Compact binary mergers (hereafter called simply mergers) 
have recently been the focus of extensive research in the ef¬ 
forts to detect gravitational waves (GW) emitted in these 
merger events (see Faber & Rasio 2012, for a recent re¬ 
view). The decreasing orbital period of the binary pulsars 
(e.g. Hulse & Taylor 1975; Taylor & Weisberg 1989; Kramer 
et al. 2006; Weisberg et al. 2010) provides so far the best 
evidence for the existence of GW. Still, as yet, GW have not 
been directly detected. While gravitational radiation should 
be produced in various different scenarios, current GW de¬ 
tectors are most sensitive to radiation emitted in mergers, 
and hence the quest of detecting GW is closely linked with 
an understanding of these events. 

A detection of electromagnetic counterparts to these 
merger events is of utmost importance (Kochanek & Piran 
1993). Such electromagnetic signals, if observed, could help 
pave the way to conhrm the hrst GW detection and would 
reveal a wealth of information about the merger process (e.g. 
Nakar & Piran 2011). Since mergers are considered as likely 
sources of short Gamma-Ray Bursts (GRBs) (Eichler et al. 
1989), these have often been considered as possible electro¬ 
magnetic counterparts to mergers. However GRBs are highly 
beamed and are only observed when their relativistic jet 


points towards us. This decreases significantly the chances 
of a coincident GW-GRB detection. Although off-axis GRBs 
cannot be observed in gamma-rays, they are followed by a 
late time isotropic radio signal, namely the orphan afterglow, 
which may possibly be observed. These orphan afterglows 
are more promising electromagnetic counterparts in this re¬ 
spect, yet their signals are expected to be weaker than the 
radio flares from dynamical ejecta from the mergers that 
we discuss here (Nakar & Piran 2011; Hotokezaka & Piran 
2015). 

In addition to GRBs, two main electromagnetic coun¬ 
terparts have been suggested. These are the “Macronova” 
(or alternatively “Kilonova”, as preferred by a few authors) 
which results from the radioactive decay of freshly synthe¬ 
sized r-process material and thus resembles a supernova (Li 
& Paczyhski 1998; Kulkarni 2005; Metzger et al. 2010; Piran 
et al. 2013; Kasen et al. 2013; Barnes & Kasen 2013; Tanaka 
& Hotokezaka 2013; Grossman et al. 2014), and radio flares 
which arise from the interactions between mergers’ outflows 
and the surrounding ISM and are analogous to Radio Su¬ 
pernovae and GRB afterglows (Nakar & Piran 2011; Piran 
et al. 2013). The latter are produced via synchrotron radia¬ 
tion of non-thermal electrons accelerated at the outflow-ISM 
shock front. Owing to the typical velocities and masses in¬ 
volved in this ejecta (~ 0.1c and ~ 10 “^Mq respectively), 
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this afterglow component should peak in the radio band on 
a timescale of a year after the merger. So far, radio flares 
have been considered only for spherically symmetric models 
(Nakar & Piran 2011; Piran et ah 2013). Our goal here is 
to characterize the effect of departure from spherical sym¬ 
metry on the resulting radio signals. This is an important 
refinement of previous results, since merger simulations sug¬ 
gests that the outflows produced by these mergers are sig¬ 
nificantly ashperical. For a related discussion on the effect of 
non-sphericity on the Macronova signal, see Grossman et al. 
(2014). 

In this work we treat only the sub/mildly relativistic 
outflow components ejected dynamically by mergers, and 
focus on the late time radio signal produced by these com¬ 
ponents. We neglect ultra-relativistic jet components that 
may produce GRBs or other short-lived signals, as well as 
non-dynamically ejected outflows such as those that may 
arise from strong neutrino driven winds (see Hotokezaka & 
Piran 2015, for a review of these components and their ex¬ 
pected radio flares). Our method follows Piran et al. (2013), 
hereafter PNR13, which assumed a spherically symmetric 
outflow, and address the impact of nontrivial geometry on 
the resulting radio signals and it’s repercussions for the de¬ 
tectability of such signals using current and future radio 
telescopes. 

The use of an approximate method is justified in light 
of the large astrophysical uncertainties regarding these sys¬ 
tems. In particular, the circum-binary density is a key un¬ 
known that will depend on the binary’s location within it’s 
host galaxy. This density is crucial in reproducing the exact 
dynamics of the system, and variations in this density will 
affect the timescale and luminosity of the radio light-curves 
dramatically. Additionally, the exact outflows expected from 
mergers are still rather ill constrained. The dynamically 
ejected outflow depends on the exact merger scenario as well 
as on unknown microphysics and in particular on the neu¬ 
tron star’s equation of state. We use the dynamical outflows 
found by Rosswog et al. (2014), yet other groups using dif¬ 
ferent numerical methods and equations of state find slightly 
different outflow characteristics (e.g Hotokezaka et al. 2013). 
Furthermore, other mass outflows arise in mergers besides 
the dynamical ejecta considered here (see e.g. Hotokezaka 
& Piran 2015). In view of these many uncertainties in the 
relevant physical parameters, we choose to use an approxi¬ 
mate analysis instead of a full fledged numerical simulation. 
As we discuss later, our approximate method maximizes the 
impact of the non-sphericity and as such it provides an up¬ 
per limit to the full effect. Our motivation is to estimate 
the detectability of these late time radio signals in com¬ 
parison with other suggested electromagnetic counterparts 
such as GRB orphan afterglows. We therefore focus here on 
evaluating only the peak timescale and peak flux of these 
signals (these will suffice for estimating the detectability), 
and we do not presume to reproduce the exact radio flare 
light-curves. We expect that the real peak time scale and 
peak flux are somewhere between those estimtated here and 
those estimated assuming spherical symmetry. 

This paper is organized as follows - we begin in § 2 
by presenting the model for the dynamics and radiation 
of a spherical outflow. This model follows and builds upon 


PNR13. In § 3 we present the piece wise-spherical approx¬ 
imation - our approximate method for treating aspherical 
outflows. We derive and discuss some general analytic re¬ 
sults that are applicable to any system in § 4. In § 5 we 
apply our analysis to various dynamic ejecta configurations 
found by Rosswog et al. (2014). We continue in § 6 by es¬ 
timating the detectability prospects of these radio flares in 
light of the results presented in § 5. Finally, we conclude in § 
7 and discuss the implications of the findings in the context 
of merger’s electromagnetic counterpart searches. 


2 SPHERICAL SYSTEMS 


Even though we are interested in non-spherical outflows in 
this work, our method is based on results for spherical sys¬ 
tems (see §3). We begin, therefore, by presenting and de¬ 
scribing some basic features of radio-flares emitted by spher¬ 
ical systems. 

The problem at hand, namely calculating the radio flare 
light curves expected from mergers’ dynamical outflows, can 
be divided into two separate tasks. The first is calculating 
the hydrodynamics of the outflow, and in particular the dy¬ 
namics of the outflow-ISM shock front. The second is calcu¬ 
lating the synchrotron flux emitted at the shock front. We 
begin by recapitulating some essential results of PNR13 (in¬ 
cluding some basic equations that are needed for our formal¬ 
ism), and later discuss some additional properties of spher¬ 
ical outflows. 


An ejected shell of matter with energy E and an initial 
velocity vq decelerates on a time scale, tcZec- 
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where n is the circum-binary ISM number density, and rrip 
the proton mass. The homologously expanding outflow is 
characterized by the energy distribution E (^ u). At earlier 
times the ejecta expands freely, whereas at later times the 
outflow follows the Sedov-Taylor self-similar solution (Se¬ 
dov 1946; Taylor 1950). At times between tdec of the fastest 
moving shell in the ejecta, and tdec of the slowest moving 
shell we describe the approximate dynamics of the contact 
discontinuity (between the shocked ejecta and shocked ISM) 
by solving: 

^nrUpR^v^ = M{R)v'^ = E v) , (2) 

o 


to obtain the radius R(t) and velocity v(t) of the contact 
discontinuity. The leftmost equality of Eq. 2 assumes spher¬ 
ical geometry and a constant ISM density. It is correct up 
to a factor of order unity, which will depend on the exact 
structure profile of the ejecta. The dynamics implied by this 
equation is governed solely by the outflow’s angular energy 
density, dEjdQ, which is just EjAir for the spherical case. In 
the following sections we generalize this for a non-spherical 
outflow confined to a solid angle Q, and we use this term 
{dE/dO) here for convenience. 


The collisionless shock wave accelerates electrons to a 
power-law distribution dN/d'y ~ This holds for 7 > 7 m, 
where 7 m is the minimal Lorentz factor of the accelerated 
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electrons, and 2 ^ ^ 3 is the power-law index which should 

be around p ~ 2.1 — 2.5 for mildly relativistic shocks, and 
around p 2.5 — 3 for Newtonian shocks. Fraction Ce and cb 
of the total energy density at the shock front are given to the 
electron and the magnetic energy respectively. Synchrotron 
radiation is emitted by these accelerated electrons. Following 
PNR13, the typical synchrotron frequency Vm and typical 
flux Fm can be written as: 

Vm{t) « 1.3 GHz (3) 




d 


10^'^cm 


^(i)> 

( 4 ) 


where d is the distance between the source and the observer 
(for simplicity we neglecting any cosmological effects). 


The synchrotron emission is self absorbed below Va , the 
synchrotron self-absorption frequency. For simplicity we fo¬ 
cus on a single standard scenario in which < i^obs, 

where i^obs is the observed frequency. This generally holds 
for jyobs = 1-4 GHz, a frequency for which numerous fa¬ 
cilities are available, as both Um and Ua will normally be 
significantly smaller than 1 GHz (see Eq. 3 and equation 5 
of PNR13 plugging in typical velocities of the merger’s dy¬ 
namic ejecta, /3 0.1 — 0.2). In this case Eqs. 3 and 4 are 

sufficient to determine the light curve. Since only the veloc¬ 
ity shells which have already decelerated at time t contribute 
to the signal, the light-curve scales as 


( 5 ) 

where the last proportionality is given by relating R{t) to 
the ejecta’s energy distribution via Eq. 2. It describes the 
fact that only the energy of shells initially moving faster 
than v(t) contribute. The complete light-curve is therefore 
given by solving Eq. 2 for the shock radius and velocity and 
substituting these into Eq. 5. 

In the following we discuss some properties of light- 
curves arising from this model. In particular we examine 
the peak flux and peak timescale for these signals which 
are important for the rest of this work, but also quantify 
the shape of the light-curve by defining the peak width. We 
begin by examining an idealized toy model for the outflow 
energy distribution E i’), namely the single-velocity shell, 
and continue by showing that any arbitrary energy distribu¬ 
tion can be approximated by this toy model for the purposes 
of evaluating the peak light-curve properties. 


A spherical “single-velocity shell” is a spherically sym¬ 
metric outflow with energy distribution dE/dv — Eq5 (I’o)- 
While this distribution does not well describe the merger’s 
outflows, we will show later on that results for single-velocity 
shells can be very useful in approximating the peak light- 
curve properties (see § 4). 


The light-curve of a single-velocity shell peaks at the 
deceleration timescale: 

\dQ,) [dfl. 
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The corresponding peak flux can similarly be expressed as: 
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The light-curve of a single-velocity shell scales as 
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Using this result we find the times at which the flux reaches 
half it’s peak value and 

A~/l = ~ 0.8 ipeafc, (9) 


ffll = « 1.5 tpeak, ( 10 ) 

where in the last equality in Eq. 10 we have substituted 
p = 2.5. With these expressions, we define the peak width. 
At, as the time period during which the flux exceeds half 
it’s peak value. Eor a spherical single-velocity shell, we find 
At = 0.73 tpeak- We use this measure for the 

peak width to quantify the extent with which two distinct 
light-curves overlap. 

An important note regarding At given above, is that 
the light-curve scaling laws on which it is based are self¬ 
similar solutions towards which (or from which) the system 
steadily evolves, and are only applicable for a single veloc¬ 
ity shell. Eor more ‘realistic’ shells which have an arbitrary 
E v)^ the light-curve around the peak flux differs signifi¬ 
cantly from this approximation, changing smoothly between 
these scaling laws which only hold at very early times or late 
times. This effect significantly increases the peak width for 
non single-velocity shells. Still, At is useful in that it may 
be used as a (rather loose) lower bound on the peak width 
of an outflow. 

We turn now to examine an arbitrary energy distribu¬ 
tion E v) outflow. The velocity Vpeak at which the light- 
curve peaks for such an outflow is easily found by differen¬ 
tiating Eq. 5 with respect to v. This peak velocity satisfies: 


d\ogE{'^ v) 
dXogv 


'^peak 


5p- 7 
2 


( 11 ) 


This result is straightforward once Eq. 5 is established - 
at earlier times than tpeak the shock velocity is larger than 
Vpeak, yet the amount of energy available E{'^ v) decreases 
faster than hence the overall flux is smaller. At 

later times, the shock velocity is smaller than Vpeak and the 
energy E(^ v) increases, yet this increase is slower than 
^-(5p-7)/2^ hence the decrease in velocity dominates and 
the total flux is once again smaller. Eigure 2 depicts a rep¬ 
resentative energy distribution for one of the merger’s dy¬ 
namic outflow configurations discussed below (nsI4nsI4), 
and shows that the velocity Vpeak at which the light-curve 
peaks corresponds to the value expected from Eq. II. 


We find that the peak flux of a radially structured out¬ 
flow can therefore be characterized using only two parame¬ 
ters Vpeak, and E(^ Vpeak)- R docs not depend on the entire 
distribution E{'^ v). Unfortunately we cannot characterize 
the time at which the light-curve peaks using only these two 
variables, for that the detailed energy distribution must be 
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Figure 1. Cumulative energy distribution E{'^ v). The red cross 
signifies the velocity at which the calculated fiux peaks. The 
dashed red line shows that the logarithmic slope at this peak 
velocity is — (5p—7)/2, as expected by Eq. 11. The cumulative dis¬ 
tribution decreases faster than this slope above Vp^aki but not as 
steeply as a step function (which would describe a single-velocity 
shell). 


accounted for. Still, we may find tpeak up to a factor of order 
unity using only these two parameters. 

We can change the time at which the Vpeak velocity 
shell decelerates while keeping Vpeak) fixed by allocat¬ 
ing more/less energy into higher velocity shells. If a large 
fraction of E{^ Vpeak) is put into very fast shells (which 
decelerate quickly), only a small fraction is left in the Vpeak 
shell, and by Eq. 1, this velocity shell decelerates earlier. An 
upper bound on this deceleration time is therefore given by 
the limiting case in which all the matter above this velocity, 
Vpeak)^ moves with Vpeak (as a single-velocity shell), 
and the distribution E{^ v) decreases infinitely fast above 
Vpeak ^ as a step function. In reality the distribution will de¬ 
crease gradually, some of the energy will go into faster shells 
with short deceleration times, and the shell moving with 
Vpeak will be affected earlier. Thus this approximation will 
be altered by a numerical factor < 1 but of order unity, that 
depends on the exact energy distribution. The peak time is 
bound from above by 

+ ^ Vpeak)\^^^ - 5/3 ^ 

tpeak ^ Uec - ( j • (12) 

A lower bound on this numerical factor (and thus on the 
timescale) may be obtained by considering the case where 
E{^ v) falls off as slowly as possible above Vpeak- Eq. 11 in 
conjunction with the fact that both E{'^ v) and dlogE/dlogv 
decrease monotonically, implies that E('^ v) must falloff 
steeper than , Hence we may find the peak time 

in the limiting case where E{^ v) decays as a power law 
above Vpeak with a slope —(5p — 7)/2. By solving Eq. 2 in 
case of a general power law profile E{^ v) oc we 

can write the expression for the shock velocity (as given by 
equation 16 of PNR13). We are interested in the case of 
k — (5p + 3)/2. The peak time is thus bound from below 
by the time at which v(t) for this power law profile equals 


Vpeak • Solving this condition by substituting the appropriate 
k and replacing Vmin by Vpeak in equation 16 of PNR13, we 
find that 


fc - 3 / 3.E(^ Vpeak) 
k \ Airrimp 


1/3 


.- 5/3 

peak 


5p-3 
5p + 3 


tdec . (13) 


Since both the upper and lower bounds scale equally 
(as tdec), we have pinpointed the peak time to within a fac¬ 
tor of {bp — 3)/{bp + 3) 0.6. We can therefore reasonably 

characterize a complicated outflow distribution by a char¬ 
acteristic velocity Vpeak (which is defined by Eq. 11), and a 
characteristic energy E{^ Vpeak), and reduce realistic radi¬ 
ally structured ejecta into roughly equivalent single-velocity 
shell building blocks. These ‘equivalent’ single-velocity shells 
(characterized by velocity Vpeak and energy E{'^ Vpeak)) will 
not reproduce identical light-curves, but will reproduce the 
correct peak flux, and approximately the correct peak time. 
Since these are the focus of this present work, such a reduc¬ 
tion will suffice when necessary. 


3 THE PIECEWISE SPHERICAL 

APPROXIMATION 

So far we have discussed spherical outflows. Turning now 
to non-spherical outflows we use a Piecewise Spherical Ap¬ 
proximation in which we do not calculate the full three- 
dimensional (3D) hydrodynamics of the system in detail. 
Instead, we decompose the 3D system into many effective 
one-dimensional (ID) radial problems which are solved sep¬ 
arately. We then combine the one dimensional results to es¬ 
timate the overall outflow evolution and the corresponding 
synchrotron emission. 

The method follows schematically the following stages: 
(i) We divide the system into solid angle elements and ex¬ 
tract the cumulative energy distribution E {'^ v) /Et for each 
element. This distribution determines the dynamics of each 
element, (ii) We calculate the spherical equivalent dynam¬ 
ics of each solid angle using E {'^ v) /Et according to Eq. 2. 
(iii) We obtain the resulting synchrotron flux for each ele¬ 
ment. (iv) We estimate the system’s overall light-curve by 
summing the contribution of all elements. 

We call this approximation “piecewise-spherical” in the 
sense that we decompose the aspherical outflow into several 
solid angle ‘pieces’ (denoted by index i) and treat each solid 
angle component i as if it is a part of a spherically sym¬ 
metric outflow with an energy distribution dE{^ v)/dEt = 
Ei{'^ v)/ELi. Using the results of § 2 for spherically symmet¬ 
ric outflows, we have a clear prescription of how to obtain 
the dynamics and synchrotron flux emitted by each com¬ 
ponent’s spherical equivalent system, and summing these 
together yields the overall light-curve under this approxi¬ 
mation. The crux of the approximation is the fact that we 
treat each solid angle component separately and indepen¬ 
dently from it’s surrounding ejecta components. This sim¬ 
plifies enormously the dynamics and the calculations but it 
has its pitfalls. In the following subsection we describe the 
limitations of this approximate method. 
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Figure 2. A cartoon sketch illustrating our method. In step (i) we divide the spherical system into N solid angle bins and 

extract the cumulative energy distribution E{'^ v) from each. In step (ii) we use each bin’s energy density distribution to calculate the 
temporal evolution of the spherical equivalent system via Eq. 2. Using these dynamics, we calculate in (iii) the flux emitted by each solid 
angle. Finally we arrive at the total light-curve by summing up the N contributions in (iv). 


3.1 Limitations of the spherical-piecewise 
approximat ion 

We turn now to discuss the important issue of the limitations 
of this approximate method. We begin by stating clearly the 
underlying assumptions of the piecewise-spherical approxi¬ 
mation: 

(1) Each solid angle element evolves according to it’s 
spherical equivalent system, and therefore the material in 
each solid angle expands only in the radial direction. 

(2) Each solid angle element evolves separately, and is de¬ 
coupled from any interactions with other elements. In partic¬ 
ular, neighboring solid angles elements do not transfer mass 
or energy from one to another and they don’t affect each 
other’s dynamics. 

Both assumptions will hold well if the angular variations 
in energy density and velocity are small, so that the system 
is nearly spherically symmetric. To illustrate this consider 
the limiting case of a jet confined within a solid angle Q. In 
this case the energy density and velocity change abruptly 
and dramatically on the boundary and assumption (1) will 
break down because the jet will expand azimuthally. This 
jet expansion occurs as the shocked region is heated and 
the hot shocked material expands into the surrounding cold 
ISM at roughly the speed of sound. If, on the other hand, 
we consider an angular region Q which is a part of a nearly 
spherical system, such a sideways expansion is hindered. 

Assumptions (1) and (2) are closely related to each 
other, as any azimuthal expansion of one bin essentially in¬ 
volves mass transfer to the adjacent bins. This mass trans¬ 
fer would in turn affect both bins’ evolution and assump¬ 
tion (2) will break down. However, we expect that even in 
case the system is far from spherical, both these assump¬ 
tions will hold reasonably well until the time of peak flux. 
This is expected because the peak flux time in this regime 
{h'm^h'a < i^ohs) IS thc deceleration time, tcZec, of the out¬ 
flow. Because tdec is by definition the time at which the 


surrounding ISM begins affecting the ejecta’s dynamics, it 
is reasonable to expect that the ISM-ejecta interactions will 
not result in any significant sideways expansion of the ejecta 
before this point. 

Faster moving ejecta components will expand az¬ 
imuthally before their slower surroundings, and energy 
will transfer from high to low velocity regions, effectively 
smoothing out the anisotropy in the initial distribution. 
Thus, our approximation overestimates the effect of as- 
phericity, and we expect the piecewise-spherical method to 
provide an upper limit, such that the true light-curve lies 
somewhere in between our present results and those calcu¬ 
lated using the spherical approximation. 


4 ANALYTIC ESTIMATES 

We begin by using the piecewise spherical approach de¬ 
scribed in § 3 to obtain simple analytic estimates that outline 
well the trends we observe later in the numerical simulations. 
To this end we compare simple asymmetric outflows with 
their spherical equivalents, defined as having the same total 
energy and total mass. In particular we examine the effect of 
asymmetry on the peak time and peak flux of the radio sig¬ 
nal. A useful picture in thinking about this problem, is that 
the transition from a spherical system into an aspherical one 
can be thought of as a “redistribution” of the system’s to¬ 
tal mass and energy into various solid angle components. 
In other words, a spherical outflow with mass Mtot and 
energy Etot uniformly distributed over a solid angle of Att 
can be morphed into an aspherical outflow by redistributing 
the mass/energy unevenly between different solid angle ele¬ 
ments. Under the piecewise spherical approach, each of these 
solid angle elements can then be evolved independently, al¬ 
lowing an estimate of the overall resulting light-curve. To 
obtain an analytic intuition to the behavior of non-spherical 
systems we consider a simple model in which we approxi¬ 
mate each solid angle element as a single-velocity shell, as 
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opposed to treating the full energy distributions E{^ v) for 
each solid angle element. 

Consider two simple examples: (a) The entire mass and 
the entire energy of the system are contained within half the 
sky (27r rad^). (b) The mass is distributed spherically but the 
entire energy is within 2tt rad^. In both cases, half the sphere 
contains no energy and therefore it does not contribute to 
the signal whatsoever. In case (a) the specific energy (and 
hence the velocity) of the relevant half sphere is the same 
as in the spherical case, yet the energy density has doubled. 
Therefore, the timescale increases by 2^^^, and the flux re¬ 
mains the same (the specific energy as well as the mass in 
half the sphere are identical to the spherical specific energy 
and total mass). On the other hand, in case (b) the relevant 
half spehere (the part which will contribute to the overall 
signal) contains the same mass density as in the spherically 
symmetric case, but twice the energy density. In this case 
the timescale decreases by 2“^^^, and the flux increases by 
2(5p-7)/4 illustrates the fact that the timescale and 

the flux of an asymmetric system may either increase or de¬ 
crease with respect to it’s spherical equivalent depending on 
the details of the mass and energy redistribution. 

We now turn to apply the same considerations to a set of 
N solid angle elements characterized by dMi , dEi , dQi where 
the index i enumerates the various solid angle elements. 
These parameters satisfy: '^■dMi = Mtot, ^idEi = Etot 
and 5]). dQi = Att. Since we are primarily interested in the 
total light curve’s peak timescale, we can estimate this under 
the piecewise spherical approximation by averaging over the 
peak flux times of all constituent solid angle elements. The 
mean peak time is a sensible estimate for the overall peak 
timescale as long as the temporal density is large enough. 
In other words, if the constituent light-curves of each solid 
angle overlap significantly around their peak fluxes, this esti¬ 
mate will be valid. This can be quantified by stating that the 
peak times of all solid angle components should be within 
~ At of each other (see Eqs. 9,10). 

We define the flux weighted mean peak time normalized 
by the original peak time of the spherical equivalent system 
as a dimensionless function r: 


T(d^^,dE,dM) = 


{tpeakj) 
Asphr) 
peak 


5p-3 _ 5p-7 

dE. 4 dM. 4 

5p-3 5p-7 

dE, 4 dM, 4 


(47r)-i/3£-J/2^y/ 


_ 5p —5 _ 


Y^dhj ^^^dEj ■* dMi 


15p-31 


EdEj 
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dMi 


where in the last equality we have introduced the normalized 
variables Xi = Xiixtot (so that 0 < < 1 and '^Xi = 1). 

This constraint eliminates three variables so that we are left 
with a total of 3 x (N — 1) free parameters. 

Particularly interesting choices of parameters are “equal 
distribution” points which retain the system’s spherical sym¬ 
metry (yielding r = 1). Equal distribution points (e.d.) ex¬ 
ist besides the trivial case of dxk = 1/^, as long as the 
mass and energy densities remain identical to the spherical 


distribution, i.e. dEkldO^k = dMkldQ^k — 1- These equal 
distribution points are critical points of r, since 


dr 

Wk 


dr 

dr 

e.d dEk 

e.d ~ dMk 


= 0 . 


(15) 


This extremum is a saddle point of r and hence with appro¬ 
priate choice of variables one can shorten or lengthen the 
timescale as desired (r < 1 and 1 < r respectively). We il¬ 
lustrate this fact graphically in Eig. 3 by plotting contours 
of T for the case of N = 2. 


We focus now on a particular scenario in which mass 
and energy are redistributed under the additional constraint 
that dE/dM is kept constant. Under this constraint, all mass 
elements have equal velocity. In this case the equal distribu¬ 
tion point discussed previously becomes a global minimum 
of T. The timescale of a constant specific energy redistribu¬ 
tion can therefore only increase with respect to the spherical 
distribution: ^ 1. 

We show this fact graphically for N = 2 in Eig. 3. 
The dashed black line in these figures depicts the curve 
dE/dM — 1. Under the constant specific energy constraint, 
T should be along this curve. It can be seen that the curve 
always passes through the equal distribution point at which 
T = 1. As one moves away from this point along the con¬ 
straining curve (in either direction), the timescale only in¬ 
creases (r > 1). 


An additional interesting point regarding the constant 
specific energy redistribution is that the total light-curve’s 
peak flux in this case should remain roughly the same as its 
spherical equivalent. Assuming the underlying time depen¬ 
dent fluxes of each solid angle overlap sufficiently around 
their peak times (under the same conditions that render 
the mean peak time a reasonable estimate of the overall 
timescale), we may use Eq. 7 to write the total light-curve’s 
peak flux as the sum of it’s components’ peak fluxes, i.e. 


Epeak ^ f dE 


5p-7 

4 

dEi = 1 . 


(16) 


Of course this is a rough order of magnitude estimate that 
can be used as an upper bound for Epeak, since the under¬ 
lying peak fluxes should not perfectly overlap. 


We have therefore shown that in case the specific en¬ 
ergy is kept constant in a redistribution of mass and energy, 
the optically thin synchrotron peak time (as estimated by 
the weighted mean r) will necessarily increase, and the peak 
flux will not change with respect to the spherical equivalent 
light-curve. Of course one can always redistribute the mass 
and energy of the system creating a fast moving jet compo¬ 
nent. Such a jet will peak at earlier times than the spherical 
equivalent system and may even lower the average peak time 
(so that T < 1). This does not contradict our results since 
this type of redistribution does not obey the constraint of 
constant specific energy. In fact, such redistribution repre¬ 
sents the exact opposite regime, in which there is a large 
(possibly extremely large) specific energy component and a 
second low specific energy component. Moreover, it is likely 
that in this case the weighted mean, r, breaks down as an es¬ 
timate of the overall timescale, as the jet component’s peak 
flux may exhibit itself as a completely separate component 
in the light-curve (in other words the two components will 
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Figure 3. The weighted mean peak time (normalized with respect to the spherical peak time) - r, as a function of the redistribution 
variables dQ, dE, dM^ for N = 2 angular divisions. Plotted are contours of r in the (dE, dM) plane, where successive figures (from top left 
to bottom right) depict decreasing values of dO.. The equal distribution point can be spotted in Fig. 3(a) at (dr2, dE, dM) = (0.5, 0.5, 0.5) 
from which it is clear that this is a saddle point. Equal distribution points exist in the following plots as well (e.g. in Fig. 3(c) the 
point is at (dQ, dE, dM) = (0.2, 0.2, 0.2)), and the value of r at these is always 1. The dashed black line corresponds to the constraint 
dE/dM = 1, which physically means that the specific energy (and hence the velocity) is kept fixed. This line always passes through the 
equal distribution point, which serves as a minimum along this line. Therefore, with the additional constraint keeping dE/dM fixed, the 
timescale necessarily increases (r ^ 1). 


not significantly overlap). As long as the dynamic ejecta re¬ 
tains a roughly spherical distribution, one should not expect 
such a jet component or any other strong variations in spe¬ 
cific energy. In this case our results stand and we expect 
a longer overall light-curve timescale. In particular we will 
show in § 5 that this is the case for the configurations we 
have studied and this serves as a motivation in developing 
these results. 


5 NUMERICAL ANALYSIS OF CALCULATED 
EJECTA DISTRIBUTIONS 

We turn now to apply our method to the results of neutron 
star merger outflow simulations by Rosswog et al. (2014), 
from which we obtain the 3D hydrodynamic configuration 


of the dynamical ejecta. These SPH simulations are a contin¬ 
uation of previous works (Rosswog 2013; Piran et al. 2013) 
which simulated neutron star mergers till several millisec¬ 
onds after the merger event, and find that roughly 10 
is dynamically ejected from the mergers, with a strong de¬ 
pendence on any asymmetry in the binary masses (equal 
mass mergers do not tidahy disrupt as much mass as non¬ 
equal binaries). Rosswog et al. (2014) followed up by simu¬ 
lating the long term evolution of the dynamical ejecta, tak¬ 
ing into account the effects of nuclear heating on the outhow 
dynamics. The major hydrodynamic effect of this nuclear 
energy deposition is in accelerating slightly the outhow and 
diminishing significant irregularities in the outhow, essen¬ 
tially smoothing out the ejecta’s angular distribution (see 
Rosswog et al. 2014, for more details). Since our calcula¬ 
tions attempt to characterize the affects of asymmetry on 
the ejecta’s afterglow signal, it is important we take the 
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most realistic initial angular distribution as conditions for 
our calculations. 

We use snapshots of the dynamical ejecta, taken from 
Grossman et ah (2014) ~ 10 s after the merger event, as ini¬ 
tial conditions for our present calculations. We choose snap¬ 
shots taken several seconds after the merger, as at these 
times the ejecta’s interaction with the surrounding ISM 
is negligible and hence the calculations of Rosswog et al. 
(2014), which did not account for these interactions, are 
valid. On the other hand, these times are late enough so 
that the nuclear decays have deposited practically all their 
energy into the outflow and affected it’s dynamics. Rosswog 
et al. (2014) show that from this point onward the ejecta 
continues expanding homologeously (neglecting ISM inter¬ 
actions), and therefore the nuclear decays have little effect 
on the outflow at later times. 

We examine outflows from simulations of three partic¬ 
ular systems, distinguished by the neutron star masses: a 
canonical equal mass 1 . 4 M 0 - 1 . 4 M 0 binary, a slightly asym¬ 
metric binary with masses 1 . 4 M 0 - 1 . 3 M 0 , and an extremely 
asymmetric system with masses I. 6 M 0 -I. 2 M 0 . We denote 
these systems ‘nsl4nsl4’, ‘nsl4nsl3’, and ‘nsl6nsl2’ respec¬ 
tively. We additionally restrict our focus in this work to the 
case where Uobs = 1.4 GHz at which numerous radio facili¬ 
ties operate. This frequency is high enough so that we are 
always in the optically thin, < b'ohs, regime. Through¬ 

out this work we take the microphysical system parameters 
as: Ce = 0.1, cb = 0.1, and p = 2.5, and use an ISM num¬ 
ber density of n = 1 cm“^, but these parameters are easily 
changed with qualitatively well known effects. 

This section is divided into two subsections loosely cor¬ 
responding to stages (i) and (ii-iv) of our method respec¬ 
tively (see Fig. 2). In § 5.1 we discuss the division of the 
systems into angular bins, and present a reduction of the ini¬ 
tial data (after binning) in terms of the characteristic energy 
and velocity of each bin. We continue in § 5.2 by presenting 
the calculated light-curves and discuss the results. Our main 
finding is an increase in the timescale of these light-curves 
in comparison to their spherical equivalents. This can be 
expected from the analytic estimates presented in § 4. 


5.1 Angular Binning 

Using spherical coordinates centered on the merger location, 
we divide each system’s outflow into iV = 20 x 20 solid an¬ 
gles, keeping the number of SPH particles in each solid angle 
fixed. In Appendix A we show that this arbitrary binning 
allocation does not affect our results. 

Figure 4 shows mollweide projections of the character¬ 
istic energy density dE (^ Vpeak) /dVt and characteristic ve¬ 
locity Vpeak of cach angular bin for the three systems. These 
are a reduction of each bin’s cumulative distribution E{^ v) 
into a single-velocity shell equivalent which reproduces the 
same peak flux and the approximate peak time (see § 2 ). 

The plots illustrate how the angular distributions tend 
toward symmetry as the progenitor merger masses become 
comparable. Still, even in the case of an equal masses merger, 
the distribution is far from spherical and there is a strong 


concentration around the merger plane (taken as ^ = 0 ) 
so that the outflow is somewhat ‘flattened’ (see also Fig. 6 
of Rosswog et al. 2014). The qualitative difference between 
equal and non-equal mass mergers is mostly apparent in 
that non-equal mergers show a stronger asymmetry inside 
the equatorial (merger) plane. This breakdown of the axial 
symmetry is a trace of the tidally disrupted tails produced 
during the merger events. As expected unequal mass binaries 
tidally disrupt each other more significantly than equal mass 
systems. 

The azimuthal symmetry breaking is especially promi¬ 
nent in the case of the nsl6nsl2 merger shown in Fig. 4(b). 
Here a large jump in the characteristic velocity is apparent 
at p 27r/3, with no coinciding jump in the energy density. 
In fact, while the (/?~27r/3, bins reach the largest 

velocities of the distribution, the same bins attain only mod¬ 
erate energy density values. The reason for this strange di¬ 
chotomous behavior is that these bins catch on to the very 
end of the tidally disrupted tail ejected by the merger. This 
tail consists of a relatively low mass fast moving outflow, 
and hence insignificant energy density. As we decrease p 
further from this point, we find a large gap region contain¬ 
ing no ejecta. This gap region is contained within the first 
(0 < 9 ^ ^ 27r/3) bin, hence the sudden drop in Vpeak- There¬ 
fore this unsmooth transition has it’s roots in the asymmetry 
of the underlying equatorial distribution. 

To illustrate this point further we plot the ns 16ns 12 
SPH particle distribution projected onto the equatorial 
plane in Fig. 5(a). The particles’ velocities at this stage are 
very nearly proportional to their radii (as they are freely 
expanding) so that the distance from the center is a tracer 
of the velocities. The angular binning in p is also plotted in 
this figure, allowing comparison with the mollweide projec¬ 
tion of Fig. 4(b). The red lines indicate the two leftmost bins 
in Fig. 4(b) over which the velocity jump occurs. This figure 
clearly shows that one bin is dominated by the end of the 
tidally disrupted tail, while the other bin is dominated by 
slower moving ejecta and encompasses the large gap region 
in the distribution (hence the significant angular size of this 
bin). 

We show that the gap region in the equatorial distri¬ 
bution does not affect the resulting light-curve significantly, 
and therefore does not introduce any binning bias. This is 
shown by excluding the gap region from the analysis and 
comparing the results with and without this region. Figure 
5(b) illustrates the SPH particle distribution in the equato¬ 
rial plane for the renewed analysis. The shaded region cor¬ 
responds to the excluded angle encompassing the gap in the 
distribution, and the radially extending grey curves show 
the renewed binning of the system. The system is rotated 
clockwise by 1.95 rad with respect to the original distribu¬ 
tion (depicted in Fig. 5(a)) so that the excluded region ends 
at (/? = 27r and the first bin begins at (/? = 0. Figure 6 (a) 
continues by plotting the characteristic velocity map of this 
system, where the black region corresponds to the gap which 
has been excluded from the analysis. The velocity in this 
case varies smoothly and monotonically from right to left 
as the bins are dominated by faster moving ejecta. Finally, 
Fig. 6 (b) depicts the light-curves for the original distribu¬ 
tion, and the distribution from which the gap region has 
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dEif^Vp^df) V erg I 
dCL V radP \ 



(a) nsl6nsl2 - characteristic energy density 


dE{f>Vpgaf) V erg' 1 

dVt L rad? J 



(c) nsl4nsl3 - characteristic energy density 


'^peak [c] 



0.1 0.12 0.14 0.16 0.18 0.2 0.22 

(b) nsl6nsl2 - characteristic velocity 


'^peak [c] 



0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26 0.28 0.3 


(d) nsl4nsl3 - characteristic velocity 





(e) ns 14ns 14 - characteristic energy density (f) ns 14ns 14 - characteristic velocity 

Figure 4. Mollweide projections of the velocity and energy density angular distributions for nsl6nl2 (a,b), nsl4nsl3 (c,d) and nsl4nsl4 
(e,f). The energy of the ejecta is concentrated in all three cases into a region close to the binary orbital plane = 0). Additionally, the 
ejecta energy is concentrated in smaller angles for larger binary mass ratios. This is a tracer of the tidally disrupted tail, which is more 
prominent with increasing binary mass ratio. While the characteristic velocity varies by a modest factor of ~ 2, the energy density varies 
by nearly 3 orders of magnitude and it is therefore the decisive factor in determining the peak timescale. 
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Figure 6. (a) A velocity map of the nsl6nsl2 system, excluding 
the gap region. The black area corresponds the excluded region. 
This distribution is rotated with respect to Fig. 4(b) (see Fig. 
5(b)). The figure illustrates that the velocity varies smoothly and 
monotonically from right to left apart from the jump across the 
excluded black region, (b) A comparison of the nsl6nsl2 light- 
curves produced by the original distribution and by the distri¬ 
bution excluding the gap region. The two light-curves are nearly 
identical, indicating that the gap region does not affect the re¬ 
sults. 




X cm 


10 

xio'' 


'^peak [c] 



Figure 5. The SPH particle distribution for nsl6nsl2 projected 
on the equatorial (x-y) plane. Ejecta velocities are proportional to 
their radius. The radially extending grey lines correspond to the 
angular binning of the system in cp. (a) The original distribution. 
The red lines emphasize the bins between which a significant ve¬ 
locity jump occurs (see Fig. 4(b)). This plot illustrates the fact 
that the velocity jump is due to the transition between the end of 
the tidally disrupted tail and a region containing slower moving 
ejecta, (b) The rotated distribution for the renewed analysis. The 
shaded black area indicates the gap region excluded from this 
analysis. 


been excluded. It is evident from this figure that the gap re¬ 
gion is inconsequential to the results, and that including it in 
the original analysis does not affect or bias the results. This 
is clear as the two light-curves presented in Fig. 6(b) overlap 
such that they are nearly indistinguishable. Additionally, we 
show more generally in Appendix A that the arbitrariness 
of any binning choice does not affect the results for any of 
the three merger scenarios. 


5.2 The Resulting Light-Curves 

After examining the angular binning and the distribution of 
characteristic variables we calculate the dynamic evolution 
of each angular bin according to it’s isotropic equivalent evo¬ 
lution, and the resulting synchrotron light-curves. The light 
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(a) nsl6nsl2 


ns14ns13 Light-Curve 



(b) nsl4nsl3 


ns14ns14 Light-Curve 



(c) nsl4nsl4 


Figure 7. The light-curves for the three binary merger systems. The flux is calculated at Fobs = 1-4 GHz (which satisfies the relation 
Frm Fa < Fobs cill bin compoueuts) and taking a typical distance of d = 10^^ cm. The dashed grey lines correspond to the spherical 
equivalent light-curves for these systems, which we compare to calculations of the present work (red lines). In all three cases the light- 
curves peak at later times than their spherical equivalents, and attain roughly the same peak fluxes. These two effects increase in their 
significance for larger binary mass ratios, which may be explained since these produce less symmetric outffows. Both these results are 
well explained theoretically (see § 4) if the characteristic velocities do not vary signiflcantly between solid angle bins (see Fig. 4 and 8). 
Quantitative results for the peak time and ffux are presented in table 1, along with these signals detectability prospects (see § 6). 


curves are compared to their spherical equivalents in Fig. 7. 
It is not surprising that the less symmetric the underlying 
system is, the greater is the differences between the results 
and their spherical equivalents. In particular the peak time 
increases dramatically in the ns 16ns 12 system when com¬ 
pared with the timescale increase of the nsl4nsl4 merger. 
The basic reason for this result is the fact that nearly all 
of the ns 16ns 12 outflow is confined to the binary plane, 
and moreover it is significantly unevenly distributed within 
this plane, as most of the energy is concentrated in a small 
tidally disrupted tail. Thus, most of the system’s energy is 
distributed over a small fraction of the sky, leading to a 
significantly larger energy density dE/dQ when compared 
with the spherical equivalent situation (in which the energy 
is distributed evenly over An rad^). As the peak time is pro¬ 
portional to {dE/dQy^^, this immediately translates into an 
increase in the light-curve’s timescale. This reasoning works 
only as long as the characteristic velocity does not vary sig¬ 


nificantly, so that the timescale is determined almost entirely 
by the energy density. 

Indeed, when looking back at figures 4(a) and 4(b) we 
notice that while the energy density fluctuates over three 
orders of magnitude, the velocity varies by a modest factor of 
two. Thus, to a zeroth order, the main effect contributing to 
the timescale delay is the compression of the ejected matter 
into a narrow solid angle which gives rise to a dramatic 
increase in the energy density dEjdO. (or equivalently to 
the mass density dM/dQ). 

In order to verify that the timescale is mainly influ¬ 
enced by the concentration of the outflow into smaller solid 
angles, we plot in Fig. 8 the peak time of each solid an¬ 
gle’s signal versus the two independent variables contribut¬ 
ing to the peak time, dM{^ Vpeak)/dQ and Vpeak- Here, we 
choose to parameterize the peak time using the mass in¬ 
stead of energy density, since the energy density correlates 
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by definition with the characteristic velocity, whereas the 
mass should be independent of it. From Eq. 12 and 13 we 
find that the peak time scales with the deceleration time, 
i.e. tpeak OC {dM{^ Vpeak)ldQ)^^^ Figure 8(a) shows 

clearly that this relationship holds quite strongly, where the 
small deviations arise from the fact that an additional degree 
of freedom (which depends on the detailed energy distribu¬ 
tion E{'^ Vpeak)) determines the exact value of tpeak- In any 
case it is evident from this figure that tpeak lies within the 
lower and upper bounds given by Eq. 13 and 12 which are 
illustrated as dashed grey lines in Fig. 8(a). 

Figures 8(b) and 8(c) show the dependence of the peak 
time on each one of the parameters independently. Also 
shown is the Pearson Correlation between the two variables. 
Clearly, tpeak depends strongly on the mass density and 
only very weakly on the velocity distribution. This supports 
our earlier conclusion that the timescale is governed by the 
ejecta’s concentration into smaller portions of the sky, and 
not by dramatic changes in the characteristic velocities. 

Finally, we compare the numerical results with the re¬ 
sults of § 4, in which we have shown that the flux weighted 
average of each solid angle’s peak time must increase for any 
redistribution in which the specific energy is kept fixed. Fig¬ 
ure 8 shows that the peak time correlates only very weakly 
with the characteristic velocity, so that the specific energy 
is inconsequential and may be treated as fixed. This jus¬ 
tifies a comparison with the fixed specific energy scenario 
discussed in § 4. Still, in order to complete the comparison 
we must ascertain that the flux weighted mean peak time, 
T, is a valid estimate of the total timescale. In § 4 we ar¬ 
gued that T is a reasonable estimate as long as the different 
contributing light-curves overlapped significantly. This was 
quantified by estimating the necessary overlap width. At, 
for single-velocity shells. Non-single-velocity shells give rise 
to larger light-curve widths, yet we can still use the single¬ 
velocity shell results as lower bounds on At in this case. 

Figure 9 depicts the flux weighted histograms of each 
bin’s peak time for the three systems , as well as the ad¬ 
ditional relevant timescales of the problem. Also illustrated 
are arrows showing the width. At, of of the underlying light- 
curves around the weighted mean. Short/dark arrows show 
the lower bounds on this width (as given by Eq. 9 and 
10 for single-velocity shells), and the longer/lighter arrows 
show the actual measured width of a typical component’s 
light-curve. The distribution of widths is centered around 

^ ^ tpeak, and ~ 1.8 X tpeak- It is clear that 

most of the components significantly overlap and hence we 
should expect the weighted mean to serve as a reasonable 
estimate of the total timescale. This is also apparent directly 
by comparing the calculated weighted mean (depicted as a 
dashed red line in these figures) with the actual total light- 
curve peak time (depicted as dashed green lines). In all three 
cases these curves are in proximity of each other when com¬ 
pared with the spherical equivalent peak time (dashed blue 
line). Another point worth mentioning is the fact that the 
actual peak time is always greater than the weighted mean 
estimate. This is most likely caused by the asymmetry of 
the light-curve around the peak (the light-curve varies at 
late/early times as /t^ respectively). As this asymme¬ 
try is unaccounted for by the weighted mean peak time, it 



(a) 




(c) 

Figure 8. The correlation of each solid angle bin’s peak time with 
its mass and velocity distributions for the nsl6nsl2 system. Fig. 
8(a) clearly shows that the peak time scales as tdec- Figures 8(b) 
and 8(c) illustrate that the main factor determining the timescale 
for this system is the mass density rather than velocity. This is 
quantified via the Pearson Correlation between the variables in 
each case. The fact the peak time is determined mainly by the 
mass density allows a comparison with the theoretical expecta¬ 
tions in case the specific energy is fixed (see § 4). The dashed grey 
lines in Fig. 8(a) show the upper and lower bounds on tpeak, cis 
given by Eq. 12 and 13. The points are all well within these limits 
with very little scatter, indi^tggfj^^^e^jjj^tXcgi^r^ef 
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-histogram 

- Weighted (tpeak) 


Total light-curve peak 
Spherical equivalent peak 
I typical At 
I lower bound on At 


Figure 9. The flux weighted distributions of peak times for each solid angle component for the three systems. The red dashed lines 
indicate the weighted mean of these components. The green dashed lines depict the actual peak time of the total light-curve. The blue 
dashed lines show the peak time in the spherical equivalent case. Arrows depict the light-curves’ width around the weighted mean, deflned 
as the times at which the flux diminishes to half it’s peak value. Dark/short arrows show the lower bounds on this width (Eq. 9 and 
10), whereas light/long arrows show the typical widths for our systems. The flgures show that the weighted mean peak time is a good 
estimate of the total light-curve’s peak time. This is consistent with the fact that the underlying components overlap signiflcantly (their 
width’s overlap over some region). 


can be expected that the actual timescale be delayed slightly 
further than this estimate would suggest. 

Analyzing ‘realistic’ merger dynamical outflows using 
our piecewise spherical approximation, we have found an 
increase in the timescale of the merger radio light-curves in 
comparison with their spherical equivalent estimates studies 
in previous works (Nakar & Piran 2011; Piran et al. 2013). 
This shifts the peak timescales to the order of ~ 10 years. In 
contrast, the peak flux of the light-curves is comparable to 
it’s spherical equivalent and it does not change signiflcantly 
when taking into account the nontrivial outflow geometry. 
Both these results are consistent with the analytic estimates 
presented in § 4 provided that the characteristic velocities 
of various solid angles are treated as roughly constant. Nu¬ 
merical values for the peak flux and peak timescales found 
are given in table 1, which also states the characteristic ve¬ 
locities and energies for the spherical equivalent cases. 


6 DETECTABILITY 

We turn now to discuss the detectability of compact binary 
mergers’ radio flares based on these results. The number of 
radio flares observable by an instrument with detection level 
Fiim in an all sky snapshot is 

Nall-sky ~ IZVAt, (17) 

where IZ is the merger event rate, V is the detectable vol¬ 
ume and At is the time that the flux is above the detection 
threshold Fum- We can improve on Eq. 17 by taking into 
account the dependence of At on the distance to the merger 
event as follows: 

/ oo 

© (F (r) - Fiim) FAt (r) = (18) 

= {if-r r M <^0 - 
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Here L is the merger’s peak luminosity, x = F(r)/Fum, 
and we assumed isotropy and neglected any cosmological 
corrections to the merger rate lZ{r) and peak flux F{r) = 
Ll{4:7Tr‘^). At (x) is the time the flux is above the detection 
limit, and it depends on the exact shape of the radio flare 
light-curve as well as on the ratio between the peak flux F 
and the detection threshold Fum- In the simple case of a 
single-velocity shell we can write down At (x) explicitly us¬ 
ing the light-curve scaling relations given by Eq. 8 (assuming 
again that < J^obs), 

At(x)= ■ (19) 

Using this result we continue in finding the relevant time- 
scale (At): 


Plugging this expression into Eq. 18 and simplifying 
the result gives an estimate for the number of merger rem¬ 
nants observable over the whole sky. This result is similar to 
that given by Nakar & Piran (2011), which estimated that 
At ~ tdec, but is larger due to a round-off error in their 
calculations. Expressing tdec and L in terms of the single¬ 
velocity outflow energy E and velocity /3c using Eq. 4 and 1, 
we find that the the number of spherically symmetric single¬ 
velocity shell radio remnants observable at 1.4GHz is: 


7vrl.4GHz 

all-sky 


123 
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( 21 ) 


We use our own canonical values in evaluating this ex¬ 
pression, which are somewhat different than those chosen by 
PNR13. We notice that the expression depends sensitively 
on the outflow velocity: Naii-sky oc which we take to be 
slightly lower than 0.2c. Naii-sky depends also on the ejecta 
energy which should be in the range of E ~ 10^° —10^^erg for 
the outflows considered. A larger value than our canonical 
10^°erg will cause an increase in Naii-sky compared with our 
canonical value. Additionally, these results are sensitive to 
the ISM density {Naii-sky oc so that binaries merging 

in less tenuous surroundings will decrease the detectability 
as well. Based on current observation of Galactic neutron 
star binaries we expect that at least a significant fraction of 
them be located in the disk of their host galaxy, where the 
typical ISM densities are ~ lcm“^ (Draine 2011). A caveat 
to this point is that it stems from the observed Galactic bi¬ 
nary neutron star population, which is both small, and con¬ 
fined within Milky-Way alone, so that other values of n are 
quite possible and this parameter is left ill-constrained. Eor 
example, Eong et al. (2014) find n ^ 5 x 10“^ —30 by evaluat¬ 
ing the 130603B short GRB afterglow. This case illustrates 
the large uncertainty in ascertaining the circum-binary den¬ 
sity using sGRB afterglow observations, and although kicks 
could relocate neutron star binaries outside their host galaxy 
where n would be quite small, there is no clear evidence to 
rule out the galactic merger scenario. 


Using these light-curves we can estimate the detectabil¬ 
ity for each merger case by evaluating Eq. 18 numerically. 
The results are given in table 1, which also states the 
peak flux and peak time found for the different merger 
light-curves. Additionally the table includes the values of 
{At) jtpeak which are the corrections to the approximation 
At ~ tpeak taking into account the specific light-curve shape. 
These calculated values (which are always of order unity) 
should be taken with a grain of salt because these light- 
curves may not necessarily reproduce the exact light-curves 
due to approximations used (which may break down after 
the ejecta begins it’s deceleration at tpeak)- In any case, a 
conservative lower bound on Naii-sky is given by replacing 
{At) I tpeak with the single velocity shell value of 0.86 (see 
Eq. 20). This substitution is straightforward since Naii-sky 
scales linearly with {At) j tpeak- 

We also provide in table 1 the outflow velocity at tpeak 
and the energy of the outflow above this velocity for the 
spherical equivalent cases. As explained in § 2, these are the 
characteristic energy and velocity which allow us to estimate 
the peak flux and time using single-velocity shell results (by 
exchanging E, v with E (^ Vpeak ), Vpeak in the relevant equa¬ 
tions). Substituting these characteristic values into Eq. 21 
for instance, yields similar results to Naii-sky calculated nu¬ 
merically. The only reason the results do not coincide exactly 
is the uncertainty in evaluating the peak time using the char¬ 
acteristic velocity and energy. This time can be found using 
these two variables only up to a factor of order unity (and 
confined between 0.6 — 1, see Eq. 12, 13). 

Our results for Naii-sky are significantly larger than 
those estimated by previous works (Nakar & Piran 2011; 
Piran et al. 2013), even for the spherical equivalent cases. 
There are several reasons for this increase. The first is a 
calculation round-off error in these papers, which leads to 
a factor of ~ 5 difference in Eq. 21 (~ 20 as estimated in 
Nakar & Piran (2011) as opposed to ~ 110 given in this 
work, taking (At) = tpeak and the canonical variables of 
Nakar Sz Piran (2011) for the comparison). The second has 
to due with the fact that the typical energies and velocities 
of the merger outflows are different than those evaluated 
Nakar & Piran (2011). Specifically, for typical energies of 2 x 
10^°erg and velocities of 0.2c we get a further increase of (2 x 
10^°/10^^)^^/®(0.2/l)^'^® 4.6 in the detectability. Lastly, 

for the case of the non-spherical light-curves calculated in 
this work, the peak time has increased by ~ 2 compared 
with the spherical equivalent cases, so that the number of 
mergers observed in an all sky snapshot doubles. 

A central issue in considering the detectability of any 
signal is it’s identification. While Naii-sky determines the 
number of signals above the detection threshold of a given 
radio facility, it is an entirely different question whether 
these signals can be identified as merger radio flares. The 
main difhculty our present results pose to the identification 
process is the increase in time-scale compared with previ¬ 
ous works (Piran et al. 2013). The time-scales we estimate 
in this work are on the order of ~ 10 yr. This long period 
will severely influence the ability to identify these signals 
as transients and therefore associate them with merger out¬ 
flows. 

In order to identify the source signal as a transient, a 
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tpeak [yr] 

Eall —sky 

(At) 

'^peak [c] 

E (^ Vyeak) [erg] 

nsl6nsl2 

0.53 

14 

9.0 X 10^ 

1.22 



Spherical Equivalent 

0.53 

6.3 

5.2 X 10^ 

1.59 

0.187 

7.04 X 10^° 

nsl4nsl3 

0.090 

10.4 

563 

1.49 



Spherical Equivalent 

0.108 

4.4 

281 

1.32 

0.174 

1.75 X 10®“ 

nsl4nsl4 

0.084 

8.7 

391 

1.38 



Spherical Equivalent 

0.093 

5.9 

287 

1.27 

0.153 

2.13 X 10^° 


Table 1. Numerical results illustrating the detectability of the systems considered. Each row designates one of the three systems and is 
divided into a top portion stating the results of the current analysis, and bottom portion stating the results for the spherical equivalent 
systems. The results are calculated at an observed frequency of = 1.4GHz, assuming a source distance ofd= lO^^cm, ISM density 
of n = lcm“^, and microphysical parameters Ce = = 0.1 and p = 2.5. The detectability is estimated by the number of observable 

sources in an all sky snapshot, N^u-sky^ assuming a merger rate of 7^ = 300Gpc“^yr“^ and a detection threshold of = O.lmJy. 

For the spherical equivalent systems only, the characteristic velocity and energy are stated in the rightmost columns. 


significant change in it’s flux must be measured. This is usu¬ 
ally achieved by returning and pointing the radio telescope 
towards the source at several different times, but on the scale 
of a decade (on which the signal changes significantly), most 
instruments will have been upgraded or replaced by more 
advanced facilities, making any recurring measurements dif¬ 
ficult to conduct and contaminated by instrumental changes. 
An alternative method of classifying a transient event is de¬ 
tecting smaller variations in the signal which would occur 
over smaller time periods. The problem with this method is 
that a typical detection should be very close to the detec¬ 
tion threshold at a low S/N ratio, and would therefore not 
be sensitive enough to small variations in the flux. Thus, 
the only way such an approach could work were if an un¬ 
likely strong signal were detected, the chances of which are 
drastically diminished. 


6.1 Comparison with GRB Afterglows 

We conclude this section by comparing the detectability of 
merger radio flares with the detectability of GRB orphan af¬ 
terglows, which are often considered as more promising elec¬ 
tromagnetic counterparts to mergers. The well constrained 
observables of GRBs are the isotropic equivalent energy of 
the burst emitted in gamma-rays and the rate of ob¬ 

served GRBs TVqrb- Orphan afterglows on the other hand 
are produced after the beamed ultra-relativistic jet has ex¬ 
panded into a roughly spherical shell, emitting isotropic ra¬ 
diation. Afterglow luminosities and rates should therefore 
depend on the total energy of the GRB and overall GRB rate 
(not the observed rate which only include bursts pointing to¬ 
wards Earth). Assuming that the GRB emission mechanism 
is very efficient, the isotropic energy of the burst should be 
comparable to the isotropic energy emitted in gamma-rays, 
Eiso ~ E^,iso (Nakar 2007). The total burst energy and over¬ 
all rate can then be calculated based on the well constrained 
Eiso and if we assume that the initial GRB jet covers 

a fraction fb of the sky: 

Etot = Eisofb ; Egrb = EcRBfb ^ • (22) 

We use these results in estimating the number of GRB 
orphan afterglows detectable at 1.4 GHz, similarly 

to the calculation for merger radio flares. The observed rate 
of GRBs is enery dependent, and can be approximated as a 
power-law over the energy range of 10^^ — lO^^erg: TVqrb ~ 
10 Gpc“^yr“^ (Eiso/lO^^erg) “ where a 0.5 — 1 (Guetta 


& Piran 2006; Nakar et al. 2006). Plugging this luminosity 
function into Eq. 21, and using Eq. 22 we find 



where we take (At) jtyeak — 1, and — 1 since the GRB 
outflow is ultra-relativistic. 


Using an expression similar to Eq. 23, Levinson et al. 
(2002) showed that although the GRB afterglow rate in¬ 
creases with smaller beaming factors, the overall detectabil¬ 
ity will decrease due to the decrease in the afterglow energy 
(see Eq. 22). Most of the variables in Eq. 23 are well con¬ 
strained, with the exception of the beaming factor fb, which 
could be on the order of = 1/20 (by definition fb < 1). 
Such small beaming factors would lower the number of ob¬ 
servable orphan afterglows in an all sky snapshot by more 
than an order of magnitude, making a blind survey detection 
implausible. 

Gomparing this result with those calculated previously 
in § 6 we see that the detectability of GRB afterglows is two 
to three orders of magnitude lower than the detectability of 
merger radio flares. This result strengthens previous claims 
by PNR13 that radio flares produced by dynamically ejected 
material are more promising electromagnetic counterparts 
to compact binary mergers than GRB afterglows. One im¬ 
portant caveat is the issue of differing timescales. While we 
have shown previously that the radio flare time-scale is of 
the order ~ 10 yr, GRB afterglows expand relativistically 
and will most likely be dominated by low energy bursts of 
~ 10^^ erg (because of the observed luminosity function) so 
that they will peak on much shorter time-scales of ~ 1 week. 
This shorter time-scale will likely make the identification of 
GRB orphan afterglows in radio surveys significantly easier. 
Still, because of the large disparity between the detectability 
of the two signals we expect that merger radio flares will be 
easier to detect, especially in bind surveys. 


7 DISCUSSION AND CONCLUSIONS 

We review here our main results and discuss their impli¬ 
cations in a broader context. In particular we address the 
question of whether the radio flares found here show promise 
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as detectable electromagnetic signals from compact binary 
mergers. We finish by suggesting further research directions 
into this and related problems. 

Throughout this work we have studied radio signals re¬ 
sulting from the deceleration of merger dynamic outflows as 
they interact with the surrounding ISM. We have expanded 
on previous works (Nakar & Piran 2011; Piran et al. 2013) by 
taking into account the asymmetry in these outflows (which 
have previously been approximated as spherically symmet¬ 
ric). We focus on observed frequency Vohs = 1-4 GHz at 
which radio facilities are most sensitive, and at which the 
outflow is optically thin. 

Although the treatment of the outflow dynamics are 
rather crudely approximated, we expect we have captured 
the ‘essence’ of the effect in the analysis - namely that 
the asymmetric outflow conhnes the system’s energy into 
a smaller portion of the sky, increasing the isotropic equiv¬ 
alent energy, which in turn increases the peak time. Note, 
however, that our approximation maximizes the effect of de¬ 
viation from spherical symmetry and as such these results 
are a sort of a limit on this effect. 

Our analytic estimates show that asphericity can either 
increase or decrease the peak timescale, depending on the 
way mass and energy are distributed. A particularly im¬ 
portant case is when specihc energy remains constant. This 
situation describes aspherical outflows expanding with the 
same velocity in each direction, such that the asymmetry is 
caused by an uneven distribution of mass in different direc¬ 
tions. More generally, this scenario reasonably approximates 
outflows in which the mass density, and not the velocity, is 
the dominant source of asymmetry. We therefore expect this 
to be the prevalent situation as long as the outflow does not 
exhibit a fast moving jet component, and in particular we 
hnd this describes well the ejecta found in merger simula¬ 
tions. In this case the analytic estimates show an increase 
in the peak timescale with respect to the spherical equiva¬ 
lent light-curve. Additionally, the peak flux should remain 
roughly the same as in the spherical equivalent case. The 
physical argumentation for these results is rather simple - 
since the outflow mass is conhned into a small solid angle, it 
must expand further radially in order to accumulate a com¬ 
parable amount of ambient mass, this takes longer (since 
the expansion velocity is unchanged). The synchrotron flux 
on the other hand is determined only by the amount (not 
distribution) of emitting mass and the velocity, which are 
both unchanged. 

Applying our method to the ejecta distributions found 
by Rosswog et al. (2014), we hnd that the light-curve’s peak 
time increases by a factor of ~ 2 compared with the spher¬ 
ical approximation, and is typically on the order of ~ 10 
years, for an ISM density of n = lcm“^. The ISM density 
is an ill constrained parameter and could be orders of mag¬ 
nitude lower if the merger occurs outside it’s host galaxy 
disk. In such cases the timescale will be even longer than es¬ 
timated above, as tpeak oc . We additionally hnd that 

the light-curve’s peak hux remains roughly the same as in 
the approximated spherical treatment. 

The implications of our hndings are important in assess¬ 
ing the detectability prospects of these signals. The longer 


timescale is a double edged sword in addressing this issue. 
On one hand, longer timescales make the identihcation more 
difficult. This is important for both all sky searches and for 
attempts to identify radio hares following identihcation of 
a GW candidate. On the other hand, the longer timescales 
mean a larger number of observable sources in an all sky 
snapshot. The combined effect will depend on the telescope 
used and on the observation strategy as well as the unknown 
background held of radio transients that might compete with 
these events. 

Support for this work was provided by an ERG ad¬ 
vanced grant “GRBs”, by the I-GORE Program of the 
Planning and Budgeting Gommittee and The Israel Science 
Eoundation grant No 1829/12 and by ISA grant 3-10417 . 
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APPENDIX A: BINNING ERRORS 

In the following, we estimate the errors associated with the 
arbitrary choice of binning, and show that they are insignif¬ 
icant. 

To estimate this, we calculate the light-curve for differ¬ 
ent divisions of solid angles keeping the number of bins fixed 
at A = 20 X 20. We enlarge and decrease randomly the size 
of each solid angle element and with it the number of SPH 
particles it contains. As we need to keep a minimal number 
of SPH particles in each bin for the analysis to work, we re¬ 
stricted the random noise so that the smallest possible solid 
angle bin would still satisfy this requirement. We addition¬ 
ally rotated the system randomly about the (p axis so as to 
begin the division at different points (essentially so that the 
first bin will not be constrained to the same location every 
time). 

Since the main variable we have studied in this work is 
the peak time of such light-curves, we extract this parameter 
and asses it’s statistical variability as a measure of the errors 
associated with the arbitrary binning process. Fig. A1 shows 
histograms for the peak timescale extracted from analyzing 
the same data using different angular discretizations. The 
standard deviations are of order ~ 1-2%, and therefore the 
binning procedure does not introduce significant systematic 
errors. 

Throughout the analysis we have made several assump¬ 
tions and simplifications that should introduce deviations 
from realistic light-curves substantially larger than a few 
percent, and therefore it is important to use caution when 
considering the numerical errors stated above. In other 
words, the stated mean peak time for these runs is only 
an approximation of the actual peak time, and the stan¬ 
dard deviations should not be taken as error bounds on this 
peak time. Nonetheless, the qualitative results such as the 
increase of the timescale due to the asymmetry have been 
established both analytically and numerically and as such 
the are robust and valid. 
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Figure Al. Distributions of the total light-curve’s peak time, 
tpeak ^ calculated using an ensemble of randomized A = 20 x 20 
angular divisions. This illustrates the variability in the peak time 
due to the arbitrary choice of solid angle bins. As this variability 
is of order ~ 2%, our results are consistent irrespective of any 
specific angular division choice. It is important to note that the 
stated standard deviations should not be taken as actual error 
bounds on our results, since the approximations of our method 
introduce significantly larger sources of errors than these (see § 
3.1). This distribution is intended only to show that the binning 
method does not introduce systematic errors. 
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